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ABSTRACT 



Context. Models of the formation, evolution and photoevaporation of circumstellar disks are an essential ingredient in many theories 
of the formation of planetary systems. The ratio of disk mass over stellar mass in the circumstellar phase of a disk is for a large 
part determined by the angular momentum of the original cloud core from which the system was formed. While full 3D or 2D 
axisymmetric hydrodynamical models of accretion onto the disk automatically treat all aspects of angular momentum, this is not so 
trivial for ID and semi-2D viscous disk models. 

Aims. Since ID and semi-2D disk models are still very useful for long-term evolutionary modelling of disks with relatively little 
numerical effort, we investigate how the 2D nature of accretion affects the formation and evolution of the disk in such models. A 
proper treatment of this problem requires a correction for the sub-Keplerian velocity at which accretion takes place. 
Methods. We develop an update of our semi-2D time-dependent disk evolution model to properly treat the effects of sub-Keplerian 
accretion. The new model also accounts for the effects of the vertical extent of the disk on the accretion streamlines from the envelope. 
Results. The disks produced with the new method are smaller than those obtained previously, but their mass is mostly unchanged. 
The new disks are a few degrees warmer in the outer parts, so they contain less solid CO. Otherwise, the results for ices are unaffected. 
The 2D treatment of the accretion results in material accreting at larger radii, so a smaller fraction comes close enough to the star for 
amorphous silicates to be thermally annealed into crystalline form. The lower crystalline abundances thus predicted correspond more 
closely to observed abundances than did earlier model predictions. We argue that thermal annealing followed by radial mixing must 
be responsible for at least part of the observed crystalline material. 
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> 1. Introduction 



With the enormous increase in the amount of high-quality ob- 
servational data of circumstellar disks in the last few years, a 
clear picture is emerging of how these objects evolve in time 
(J0rgensen et al. 2007, 2009:'Lo onev et al.ll2007HLommen et al.1 
i2008; Sicilia -Aguilar et al. 20091) . They form during the col- 
lapse of a pre-stellar cloud core, undergo a number of ac- 
cretion events (FU Orionis and EX Lupi outbursts), live 
for 3 to 10 Myr, and, shortly before they are destroyed, 
open up huge gaps visible in the dust con t inuum and some- 
times also in gas lines (iD'Alessio et al.l 2005t Goto et all 
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2006HBrittain et all2007HRatzka et al.ll2007tlBrown et al.l2008t 



' iPontoppidan et all l2008l) . These physical changes are echoed 
in the evolution of their chemical composition and dust prop- 
erties. Pre-stellar cores contain mostly simple hydrides, rad- 
icals and other s mall molecules, largely frozen ou t onto the 
cold dust grains (iBergin & Langeil 119971: iLee et al.1 12004 . A 
fully formed disk is predicted to contain a much richer chem- 
ical mixture with a w ide variety of complex org anic molecules 
dRodgers & CharnlevI 120031: lAikawa et al.ll2008h. although only 
simple organics have been o bserved so far dLahuis et al.ll2006t 
ICarr & Naiitall2008l;ISalvk et al. 2008). The dust by this time has 
grown from less than a micron to millimetres and centimetres, 
and part of it has evolved from an amorphous to a crystalline 
structure (Bouwmanet al. 2001, 2008; van Boekel et al. 2005; 
iNatta et al..2007;,Lommen et al..2007M2009;,Watson et al..2009u 



lOlofsson et al.ll2009 *). Crystalline silicate dust is observed down 
to temperatures of 100 K, well below the threshold of 800 K re- 
quired to convert amorphous silicates into crystalline form. One 
of the central questions of this paper is how much silicate ma- 
terial comes close enough to the star to be crystallised. We also 
investigate how the crystalline silicates end up so far outside of 
the hot inner disk where they appear to be formed. 

One way to answer these questions is to construct detailed 
models of the evolution of circumstellar disks based on our cur- 
rent understanding of the physics of these objects, and then com- 
pare to the available observational data. However, it would re- 
quire extraordinarily heavy computations to run a model that 
does justice to all physical processes known to be involved. 
A circumstellar disk ranges from a few stellar radii to hun- 
dreds of AU, and lasts for several million years. A full model 
would therefore have to resolve hundreds of millions of inner or- 
bits, and span some five orders of magnitude on a spatial scale. 
Moreover, an accurate radiative transfer method is required to 
properly compute the temperatures. All this is clearly too de- 
manding. Most multidimensional hydrodynamical simulations 
therefore solve sub-problems that only capture part of the disk, 
or only evolve over a limited time. Even these models, though, 
require days or weeks of CPU time for a single set of parameters. 

An alternative approach is to parameterise most of the 
physics in some form, and treat the disk evolution as a sim- 
pler one-dimensional (ID) time-dependent problem. One as- 
sumes axisymmetry and integrates the density vertically to ob- 
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tain the surface density E, which is now only a function of 
the radial coordinate R and the time t. These kinds of models 
go back to the pioneering work b y Shakura & Sunvaev ( 1973) 
and iLvnden-Bell & Pringld (Il974l) . In order to use these mod- 
els throughout the disk's lifetime, some way must be found to 
also include the birth phase of the disk in a reasonably realistic 
way. Two-dimensional axisymmetric hydrodynamical models of 
disk formation show the presence of a stand-off shock that de- 
celerates the supersonically inf ai ling matter as it approaches the 
disk from ab ove and below (e.g., Tscharnutenll987l : lYorke et aP 
[1993 ; Neufel d & Hollenbachll 19941) . The structure of this stand- 
off shock is clearly multidimensional in the outer regions, but 
may be approximated in a simp ler manner in the inner regions 
dNakamoto & Nakagawal 11994 . iHueso & GuiUo^ ( l2005h con- 
structed a ID disk evolution model with such an approximation 
and used it to analyse two T Tauri stars. This showed that sim- 
ple models of disk formation and evolution can be very powerful 
and yield valuable insight into the evolutionary stage of young 
stellar objects. Similar models have also been used to analyse the 
statistics of the accretion rates measured in pre-main-sequence 
stars Imilemond et al. 2006b; Vorobyov & Basu 2008). 

Another problem addressed with these ID parameterised 
models is the origin, evolution and transp ort of gas and d ust in 
circumstellar disks. For instance, Dullemo nd et alj (l2006al here- 
after DAW06) suggested that the initial outward expansion of the 
disk during the disk formation phase (observationally the Class 
O/I phase) may be very effective in transporting thermally pro- 
cessed dust to the outer parts of the disk. Based on that work, 
one would expect to find a number of disks with nearly 100% 
crystalline d ust. However, n o such extremely crystalline disks 
are observed (JBouwman et al. 2008: Watson et al. 2009) . This is 
one of the issues addressed in this paper. 

Yet another application of parameterised disk evolution mod- 
els was shown in Visser et al. (2009, hereafter V09). We fol- 
lowed the envelope material from thousands of AU inwards, 
through the accretion shock and into the disk, to analyse when 
ices evaporate from and recondense onto the grains. Since much 
of the interesting physics happens in the outer regions of the disk 
(several hundred AU), where the accretion shock no longer has 
a simple ID shape, some way had to be found to include the 2D 
axisymmetric nature of this region without having to resort to 
a full-scale multidimensional hydrodynamical simulation. Our 
recipe was to construct a semi-2D disk model, i.e., to generate 
the 2D density structure p(R, z, f) (where z is the vertical coor- 
dinate away from the midplane) out of the ID surface density 
"LiR, t) using the disk temperature to compute the scale height 
at every radius. The accretion onto the disk was then modelled 
by following infalling matter on ballistic supersonic trajecto- 
ries until its density equals the density of the disk. This yields 
a 2D axisymmetric shape of the stand-off shock that is similar to 
that obtained by full 2D axisymmetric hydrodynamical models 
(lYorke et al.i ri993; Brinch et al. 2008). However, it is not triv- 
ial to link this type of multidimensional infall structure back to 
the ID disk evolution model. The main obstacle is that a sub- 
stantial fraction of the matter does not really fall onto the disk's 
surface but onto the disk's outer edge, i.e., it accretes from the 
side instead of from the top. Moreover, this matter rotates with 
a very sub-Keplerian velocity. If this is not treated in some way, 
the total angular momentum balance of the system is violated, 
potentially leading to very wrong estimates of the evolution of 
the disk mass with time. We devised an ad-hoc solution to this 
problem in V09, which gave reasonable results but did not prop- 
erly conserve angular momentum. We derive a more rigorous 
solution in the current paper. 



The problem of sub-Keplerian accretion onto a disk is not 
studied here for the first time. It has long been known that mate- 
rial falling in an elliptic orbit onto the surface of the disk has sub- 
Keplerian angular momentum, even in the case of infinitely flat 
disks. When mixing with the disk material, which is in Keplerian 
rotation, a tor que is exerted that pu s hes th e disk material in to- 
wards the star. ICassen & MoosmanI (Il98lh described this prob- 
lem and solved it elegantly. For disks without internal torque 
and with small vertical extent, they found an analytical solution 
for the disk surface density as a function of ti me, and they also 
presen ted numerical results for viscous disks. iHueso & GuillotI 
(2005) derived an alternative solution: they simply calculated 
the radius for which the specific angular momentum of the in- 
falling matter is Keplerian, and inserted the matter into the disk 
at that radius. This is not unreasonable, because one may as- 
sume that the material, after it hits the surface of the disk, first 
adjusts its own orbit before it mixes with the disk material be- 
low. It may be argued that as long as the angular momentum 
budget is not violated, the disk evolution is not much affected 
by the precise treatment. However, if one wishes to follow the 
radial motion and mixing of certain chemical species or types of 
dust, this may depend more sensitively on the way the angular 
momentum problem is treated. 

In this paper we add a consistent treatment of sub-Keplerian 
2D accretion to a semi-2D model for the formation and viscous 
evolution of a disk (Sect.|2]i. Some basic disk properties arising 
from the n ew treatment are discussed in Sect. [3] and the results 
from |V09l are b riefly revi sited in Sect. ID The dust crystallisa- 
tion results from lDAWOd are re-evaluated in Sect.|5] and finally 
conclusions are drawn in Sect. |6l 

2. Equations 

The model describes the formation and evolution of a circum- 
stellar disk inside a collapsing cloud core. The disk's surface 
density E is governed by two processes: accretion of material 
from the envelope and viscosity-driven diffusion. This requires 
the continuity equation 



dim) djYRuR) 



dt 



dR 



(1) 



with ur being the radial velocity and S the source function that 
accounts for accretion from the envelope. In addition to conser- 
vation of mass, ensured by Eq. ([T), there must be conservation 
of angular momentum ( Lvnden-Bell & Pringle 1974) : 
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with V being the viscosity coefficient (discussed further down) 
and Qk the Keplerian rotation rate. 

In the last term in Eq. (|2]i, belonging to the accreting mate- 
rial, the rotation rate Q is smaller than Qk. If accretion would 
occur with the Keplerian velocity, Eq. (|2]) can be solved to give 
the expression for ur used by DAWQ^ 



Ur(R, t) 



d 



Y^/R^R 



{I.v^/R) . 



(3) 



Working from that solution, two methods have been used in the 
recent literature to correct for the sub-Keplerian velocity of the 
accreting material. Hueso & Guillot (2005) and DAW06 modi- 
fied the source function so that all incoming material accretes at 
the exact radius where its angular momentum equals that of the 
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disk. One might picture this as a quick redistribution of material 
in the top layers of the disk after accretion initially occurred in 
a sub-Keplerian manner A disadvantage of this method is that it 
introduces a discontinuity in the infall trajectories: upon accre- 
tion onto the dis k, ma terial instantaneously jumps to a smaller 
radius. Hence, in lV09l we chose to modify the expression for the 
radial velocity instead of that for the source function, taking 






(4) 



with M, being the stellar mass at time t and 77, a parameter with a 
constant value of 0.002. While this method has the desired effect 
of transporting material to smaller radii without discontinuities, 
it does not properly conserve angular momentum. 

As an alternative to these two methods, Eq. (O can also be 
solved more generally for the case that Q < Q.^^, giving 
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Expressi ons for Q (or, rather, for the azim uthal velocity) mayb e 
found in ICassen & MoosmanI (Il98lh a nd iTerebev et all (11984 . 
The source function is taken from lV09l : 



S(R,t)^2pu,., 



(6) 



with p being the density at the disk-envelope boundary and m, the 
vertical component of the envelope velocity field at that point. 
Equations ([T) and ^ now describe a fully continuous solution to 
the evolution of the surface density with proper conservation of 
angular momentum. A physical interpretation of Eq. Q is pro- 
vided in Sect. [3] Note that the model contains no prescription for 
photoevaporation or other disk dispersal mechanisms that come 
into play in the T Tauri or Herbig Ae /Be phase. This could eas- 
ily be added (e.g., Alexander et al. 2006), but we choose to focus 
only on the infall phase of star formation. Hence, our simulations 
are only reliable up to an age of about 1 Myr 

The boundary between the disk and the envelope was com- 
puted in[V09 as the surface where the density from both was the 
same. Here, we take instead the surface where the ram pressure 
of the infalling gas equals the thermal gas pressure from the disk: 



2 _ k-l sPdisk 
Penv^ — 



(7) 



with u being the velocity of the infalling material and yu the mean 
molecular mass of 2.3 proton masses per gas particle. The tem- 
perature at the surface of the disk, T^, is determined by irradia- 
tion from the star and by viscous heating. The use of the isobaric 
surface as the disk-envelope boundary instead of the isopycnic 
surface can have some effect on the model results, but, as shown 
in the next sections, the effects of the new treatment of the sub- 
Keplerian accretion are usually much larger. 

The viscosity coefficient a ppearing in Eqs^(l2]i and (|5]l is cal- 
culated in the same way as in DAW06 and V09, i.e., with the a 
prescription from S hakura & Sunyaevi ( il973i) : 



v{R,f)^ 






(8) 



with I'm being the midplane temperature. Effects of self-gravity 
of the disk are ignored in Eq. ^, likely leading to an over- 
estimate of the disk-to-star mass ratio for the more massive 
disks (Mj <: 0.3M ,) covered by our choice of initial conditions 
(IVorobvovll2010h . For the current work, the a parameter is set 



to a constant value of 0.01 (iHartmann et alJ fl998h . Choosing 
a different value, or allowing a to vary in time or with radius, 
would affect the evoluti on of the disk (e.g., L i n & Pringle 199^ 
Kratteretal. 2008; Vo robvov & Basd l2009t iRice & Armitagel 
,2009; ,Rice et al..,2010.) . In order to keep the current work fo- 
cused on the new treatment of the sub-Keplerian accretion, the 
effects of different viscosity prescriptions will be explored sepa- 
rately in a future publication. 

3. Size, mass and structure of the disk 

The net effect of the rightmost terms in Eqs. © and (|5]l is the 
same: to provide an additional inward flux of matter. However, 
the radial dependence of the additional flux is different. In the 
old method (Eq. (HJi), it was largest at the inner edge of the disk 
because of the 7? '^- proportionality. In the new method, on the 
other hand, it is strongest at the outer edge of the disk because 
of the sharp decrease in surface density at that point. Physically, 
the new dependence is easy to understand. Near the disk's in- 
ner edge, accretion occurs into a large column of material (large 
S), and the torque arising from the different azimuthal velocities 
only results in a weak inward push. The same torque provides 
a much stronger push at the outer edge, where there is less disk 
material per unit surface (small S). In effect, the material falling 
onto the disk's outer edge pushes the disk inwards and hmits its 
radial growth. 

In order to quantify the size of the disk, we first have to de- 
fine the outer edge R^. Since the equations allow an infinitesi- 
mal part of the disk to spread to an infinitely large distance, we 
cannot simply use the full radial extent of the disk as its outer 
edge. Instead, we define R^ as the radius that contains 90% of 
the disk's mass. In general, the surface density at that point re- 
mains well above 0.01 g cm"^ (see, e.g.. Fig. |2]i, so we do not 
have to worry about the outer parts being rapidly photoevapo- 
rated or dispersed for the duration of our simulations. 

With the new method, the disk's outer edge initially lies be- 
yond the centrifugal radius. 



RM = —c.mlPSll , 



(9) 



where Cs is the so und speed, mo is a numerical factor equal to 
0.975 (Shul ll977h . and Qq is the solid-body rotation rate. The 
outer edge grows more or less linearly in time as the disk spreads 
out to conserve the angular momentum of the entire system. The 
centrifugal radius grows as t^, so the ratio Rd/Rc decreases as 
the collapse goes on. If the rotation rate is high enough, the cen- 
trifugal radius overtakes the outer edge of the disk before the 
end of the collapse. For the remainder of the collapse phase, R^ 
equals R^, because the disk cannot be smaller than the centrifu- 
gal radius ("C assen & Moosmanlll98ll) . If the rotation rate is low 
enough that the centrifugal radius does not overtake the outer 
edge of the disk before the end of the collapse, /?d continues its 
near-linear growth until face- After accretion from the envelope 
has ceased and the inward push from the accreting material is 
no longer present, the disk is free to expand more rapidly. Its 
outer edge again grows linearly in time in both the high- and 
low -rotation cases, and does so at a higher rate than in the initial 
linear-growth regime. 

Figure [T] shows t he ou ter edge as a function of time for the 
standard model from lV09l which has an initial core mass Mo of 
1 .0 Mq, a sound speed c^ of 0.26 km s"' , and a low rotation rate 
Qo of 10"''* s"'. The disk grows to about 40 AU at the end of 
the collapse phase at face = 2.5 x 10^ yr, and it spreads to ten 
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Fig. 1. Oute r disk radius as a function of time for the standard 
model from IV09i Solid black: method from the current paper; 
grey: method from ' V09'. The dash-dotted curve shows the cen- 
trifugal radius, and the dotted vertical line indicates the end of 
the envelope accretion phase. Note that R^ is a physical quantity 
only up to face; for larger f, we simply plot the same mathematical 
function (Eq. ©). 



Fig. 2. Surface density of the disk for the standard model from 
V09 at 0.15, 0.25 and 0.50 Myr (0.6, 1.0 and 2.0 face)- The grey 
lines show the results from V09; the black lines show our new 
results. The filled circles indicate the radius that contains 90% of 
the disk's m ass. The thin do tted line is the minimum-mass solar 
nebula from|Hayash| (Il981h . 



times that size over the next 7.5 x 10^ yr. The method from lV09l 
(labelled "old" in Fig. [TJ provides a more rapid growth during 
the collapse phase, giving an R^ of 230 AU at face- The post-infall 
spreading occurs at the same rate as in the new method. We find 
the same qualitative differences and similarities between t he old 
and the new method for the rest of the parameter grid from lV09l 

Is the smaller disk from the new method a realistic scenario? 
Our model necessarily assumes ur to be the same at all heights 
above the midplane. In reality, one might expect the infalling 
material to interact mostly with the surface of the disk, exert- 
ing its torque on only a fraction of the total surface density. This 
pushes the surface material radially inwards, while material near 
the mid plane is unaffected. Indeed, the hydrodynamical simula- 
tions of lBrinch et alJ (l2008h show such behaviour, with material 
near the surface moving inwards and material at the midplane 
moving outwards. However, as this midplane material moves out 
beyond the outer edge, it becomes surface material itself, and is 
in turn exposed to the inward push from the envelope material. 
The assumed altitude-independence of mr may therefore lead to 
the outer radius being underestima ted by a factor of two or t hree. 
On the other han d, the models of iRice & Armitagg (l2009h and 
iRice et aP (|2010') suggest that disks with a spatially dependent 
viscous a expand less slowly than do disks with a constant a, in 
which case we would be overestimating the outer radius of the 
disk. Altogether, the new method produces disk sizes that are 
consistent with observations and other simulations, and it offers 
an improvement over the unrealistically large disks sometimes 
produced with the old method. 

Although the disks are smaller with the new method, they 
are not necessarily less massive. In fact, we find a rather large 
mass increase for the standard model from V09: from 0.05 
to 0.13 M© at the end of the collapse phase. There are three 
causes for this, each of which accounts for 0.02-0.03 M©: (1) 
the new definition of the disk-envelope boundary (Eq. O); (2) 
the new sub-Keplerian coiTection (Eq. (|5])); and (3) an improved 
integration scheme in the computational code. The effects are 
smaller for more rapidly rotating clouds. For example, the disk 
mass obtained for the reference model from iV09i (Mq = 1.0 



Mo, Cs = 0.26 km s"', Qq = 10"'^ s"') is unchanged at 0.43 
Mq. Observed disk masses are t ypically between 0.0()1 and 0.1 
Mq for low-mass protostars (e.g. [Andrews & Williamsl2007al[bl) . 
For realistic initial conditions (i.e., Qq = 10"'"^ s""' rather than 
10"'^ s"'), our model is known to overpredict the disk-to-star 
mass ratio by about a factor of three compared to observations 
(|j0rgensen et al.ll2009h . However, it has also been argued that 
observed disk masses have been systematically underestimated 
(Hartmann et al. 2006), a view suppo rted by the hydrody namical 
simulations of iVorobvovl (120091) and lKratter et all (1201 Oh . 

The combination of smaller sizes and equal or larger masses 
means that the surface density profiles have also changed from 
the old to the new method. Figure|2]shows the surface density for 
the standard model from V09 at three different time steps. The 
time evolution is qualitatively the same for the old (grey) and the 
new method (black). Up to the end of the accretion phase (t - 
face; dashed curves), mass is mostly added onto the outer parts 
of the disk, and the surface density inside of /?d remains nearly 
constant. At later times, when there is no more accretion from 
the envelope, the disk continues to spread and the surface density 
inside of /?d gradually decreases. Quantitatively, the rate at which 
the surface density changes is different for the old and the new 
method, and also depends on the choice of initial conditions. 

At face, the surface density from our new method falls off 
as /? ' '' inside of 2 AU. This is the same slope as in the 
original model for th e minimum-mass solar nebula (MMSN; 
IWeidenschillingI [19771; lHayashl.1981.) . but it is much shallowe r 
than the slope of -2.2 in the new MMSN model of lDeschl(l2007h . 
The slope of - 1 .5 is also commonly found in other disk evolu- 
tion simulations (e.g., Vorobyov & Basu 2009; Rice et al. 2010). 
However, whereas those simulations tend to have a constant 
slope out to at least a few tens of AU, our surface density flat- 
tens off to a slope of -0.7 between 2 and 30 AU. This is due to 
the mass being accreted mostly onto the outer parts of the disk. 
Indeed, the old method, where the infalling material was dis- 
tributed more evenly across the entire disk, shows a somewhat 
steeper slope of about -1 .0 between 2 and 30 AU. 
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4. Gas-ice ratios 

During the collapse of the pre-stellar core to form a protostar and 
circumstellar disk, large changes occur in both density and tem- 
perature. Many molecular spe cies are frozen out onto dust grains 
befor e the onset of collapse dBergin & Langerlll997t iLee et all 
l2004i) . The warm-up phase during the collapse causes some of 
them to evaporate, and they may freeze out again once mate- 
rial settles near the disk's relatively cold midplane. In V09, we 
modelled these processes for carbon monoxide (CO) and water 
(H2O). Here, we investigate whether our new sub-Keplerian ac- 
cretion correction affects these results. 

The model consists of several steps. First, it computes the 2D 
axisymmetric density and velocity structure at regular intervals 
from the onset of collapse (t - 0) to the point where the en- 
tire envelope has accreted onto the star and disk (t = face)- Dust 
temperatures are computed at the same time intervals with the 
radiative ti-ansfer code RADMC (Dullemond & Dominik 2004), 
and the gas temperature is assumed equal to the dust tempera- 
ture. Using the velocity structure, the model then computes in- 
fall trajectories from a grid of initial positions in the envelope. 
Each trajectory represents an individual parcel of gas and dust, 
for which we now know the density and temperature as a func- 
tion of time and position. This allows us to compute the adsorp- 
tion and desorption rates of CO and H2O, and solve for their gas 
and ice abundances in a Lagrangian frame. Finally, the parcels' 
abundances are transformed back into 2D axisymmetric abun- 
dance profiles for the disk and remnant envelope. 

As discusse d in S ect. [3] the main difference between the disk 
properties from vOy and the new method is the size of the disk. 
Because the mass has increased or stayed the same (depending 
on the model parameters), the density of the disk is now higher. 
In the absence of viscous and shock heating (see below), this 
leads to a more rapid decrease of the dust temperature along the 
midplane at smal l radii. For example, at f = face in the reference 
model fromlVOl (Mq = 1.0 Mq, c, = 0.26 km s-\ fio = 10"'^ 
s"'), the temperature at 30 AU is 45 K. With the new method, the 
temperature is down to 3 1 K at that point. The midplane temper- 
ature decreases further with radius until it reaches 21 K at 160 
AU with the new method. At larger radii, photons scattering off 
the surface of the disk begin to reach the midplane again and 
the temperature gradually in crease s to 25 K at the outer edge at 
500 AU. The original disk in IVOgl extended to 1500 AU, and the 
turnover from a decreasing to an increasing temperature along 
the midplane occurred at 1000 AU instead of at 160 AU. Hence, 
even though the midplane temperature decreased less rapidly in 
IV09l it still reached a lower minimum: 15 K, as opposed to 21 
K with the new method. 

If all CO is taken to desorb at 1 8 K, as it would for a pure CO 
ice, the new temperature profile does not allow for any solid CO 
to exist in this particular model at face- At later times, when the 
disk has spread to larger sizes and the protostar has become less 
luminous ( V09), a region appears around the midplane where the 
temperature does go below 18 K and CO freezes out again. In re- 
ality, solid CO forms a mixture with solid H2O, and some of the 
C O remains t rapped in the ice matrix a t temperatures above 18 
K (IColUngs e t al. 2004: IViti et"aDl2004l: Fayolle et al. in prep.). 
When the gas-ice ratios are computed accordingly, the mass frac- 
tion of solid CO averaged over the entire disk at face goes from 
33% with the old method to 20% with the new method. 

T he gas-ice ratios for pure CO in the standard model from 
rV09l (Mo = 1.0 Mo, c, = 0.26 km s"', Qq = lO"''* s"') are 
unchanged. With both the old and the new method, the entire 
disk is warmer than 18 K at face, so there is no solid CO. If we 



allow part of the CO to be trapped in the H2O ice, the disk- 
averaged solid fraction is 15% with both methods. 

Massive disks (Mj > 0.3M,) are susceptible to heating 
sources other than the stellar radiation field, such as viscosity 
and gravitational shocks. These heating mechanisms are not in- 
cluded in the chemical part of our computational code, where 
the gas-ice ratios are computed. (Viscous heating is, however, 
included in the part of the code that calculates the evolution of 
the mass and surface density of the disk and the evolution of 
the crystalline silicate abundances.) Viscous and shock heating 
would result in higher temperatures for the more compact disks 
produced by our new method than they would for the more ex- 
tended disks produced by our old method. However, this is un- 
likely to change the solid fraction of CO by much. At the end 
of the collapse phase, the entire disk is already warmer than the 
evaporation temperature of pure CO (18 K), so all CO is pre- 
dicted to be in the gas phase. Adding additional sources of heat 
will only strengthen that conclusion. For CO ice mixed with H2O 
ice, the fraction of trapped CO hardly changes between the re- 
spective evaporation temperatures of CO and H2O (Fayolle et al. 
in prep.), so the additional heating from viscosity and shocks 
would have to bring the temperature to more than 100 K to 
change the CO ice abundance significantly. This will only hap- 
pen in a small zone around the snowline, so the mass fraction of 
solid CO averaged over the entire disk will remain close to the 
value of 20% mentioned above. 

In summary, the new treatment of the sub-Keplerian accre- 
tion results in disks that are a few degrees colder in the inner 
parts and a few degrees warmer in the outer parts. Overall, the 
new CO ice abundances are between and 50% smaller than 
those obtained with the old method. In all cases, H2O remains 
solid in the entire disk except for the inner few AU. 



5. Crystalline silicates 

5. 1 . Observations and previous model results 

Infrared spectroscopic observations have shown that 1-30% of 
the silicate dust in the disks around Herbig Ae/Be and clas- 
sical T Tau r i stars occurs in crysta lline form (JBouwman et al.l 
200 il I2OO8I: fvan Boekel et"ai]|2005h . Even larger fractions (up 
to 100%) are found in the inner 1 AU of some sources 
(Watson et al. 2009) . The interstellar medium, from which this 
dust originates, has a cry stalUne fraction of at most 1-2% 
dKemper et alJ l2004 l2005h . Two mechanisms are thought to 
dominate the conversion of amorphous silicates into crystalline 
form during the formation and evolution of a circumstellar disk. 
At temperatures above ~1200 K, the original grains evaporate 
(iPetaev & Wood 120051) . When the gas cools down again, the 
silicates recondense in crystalline form dDavis & Richteijr2003l : 
Gail 2004). Alternatively, amorphous dust can be thermally 
anneale d into crysta lline dust at temperatures above ~800 K 
(.Wooden et alj 120051) . As the disk is formed out of the parent 
envelope, part of the infalling material accretes close enough 
to the star that it is heated to more than 800 or 1200 K and 
can be crystallised. However, the observations show significant 
fractions of crystalline dust at least out to radii corresponding 
to a temperature of 100 K. Crystalline dust is also found in 
comets, which are formed in regions m uch colder than 800 K 
(IWooden et alJll999l: iKeller et alJl200 6). This suggests an effi- 
cient radial mixing mechanism to transport crystall ine material 
from the hot inner disk to the colder outer parts dNuthl [19991 : 
iBockelee-Morvan et alJl2002HKeller & C^l2004 . 
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An argument against large-scale radial mixing was recently 
provided by spatially resolved observations with the Spitzer 
Space Telescope. Bouwmanet al. (2008) found a clear radial de- 
pendence in the relative abundances of forsterite and enstatite, 
two common specific forms of crystalline silicate. If both are 
formed in the hot inner disk and then transported outwards, one 
would expect the same relative abundances throughout the entire 
disk. Hence, the observed radial dependence argues in favour of 
a localised crystallisation mechanism such as heating by shock 
waves triggered by gra vitational instabilities (JHarker & DeschI 
l2002tlDesch et al.l2005h . However, the observations do not com- 
pletely rule out the possibility of crystallisation in the hot inner 
disk followed by radial mixing; at best, they provide an upper 
limit to how much crystalline dust can be formed that way. 

In another set of Spitzer observations, crystalline spectro- 
scopic features in the 20-30 fim region were detected three 
time s more frequently t han the crystalline feature at 11.3 
/im dOlofsson et al.l 120091) . This is unexpected, because shorter 
wavelengths trace warmer material at shorter distances from the 
protostar, where all models predict the crystalline fractions to be 
larger The crystalline 11.3 //m feature may be partially shielded 
by the amorphous 10 //m feature, but Olofsson et al. showed that 
this alone cannot explain the observations. A full compositional 
analysis (Olofsson et al. subm.) is required to shed more light on 
this "crystallinity paradox". 

In the model of DAWOg, crystallisation occurs right from 
the time when the disk is first formed. Indeed, because the disk is 
very small at that time, its dust is hot and nearly fully crystalline. 
As the collapse proceeds and the disk's outer radius grows, an 
ever larger fraction of the infalling material does not come close 
enough to the star anymore to be heated above 800 K. In the 
absence of strong shocks, this results in amorphous dust being 
mixed in with the crystalline material. Hence, the crystalline 
fraction averaged over the entire disk is expected to decrease 
with time. There is tentativ e observational support for an age- 
crysta llinity anticorrelation (Ivan Boekel et al. ]2005 : Apai et al J 



2005h. but this is far from conclusive dBouwman et al.l l2008t 



Watson et al.ll2009 l). One should of course consider the fact that 



obse rvations d o not probe the entire disk. If the model results 
from lDAW06l are interpreted over a limited part of the disk, such 
as the 10-20 AU region, they show only a small difference in the 
crystalline fractions at 1 and 3 Myr. Add to that the uncertainties 
in the ages for individual objects, and it is clear that the model 
results cannot be said to conflict the observational data. 

The crystalline fractions obtained bv lDAW()6l were all on the 
high end of the observed range of 1-30%, unless unreasonably 
high initial rotation rates were adopted for the envelope or the 
disk temperature was lowered artificially. If we accept that only 
part of the silicates in the outer disk originate in the hot inner 
region - so the other part, formed in situ, can account for the 
ob served rad ial abundance variations - the discrepancy between 
the lDAWQq model and the observations becomes even larger. In 
the following, we show that we obtain more realistic crystalline 
fractions with our new method. 

5.2. New model results 



The two main differences between the old method of iDAWOa 
and our new method are (1) the treatment of the disk as a multi- 
dimensional object instead of just a flat accretion surface and (2) 
the improved solution to the problem of sub-Keplerian accretion 
(Sect. |2]l. The former has the largest impact on the crystallisa- 
tion. In case of a fully flat disk, material falls in along ballistic 
trajectories until it hits the midplane at or inside the centrifugal 
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lace) for the default model 

of|] 

of Re (2.2 AU) and /?ann (3.6 AU). Top: mass-loading as a func- 
tion of radius. Bottom: infall trajectories (solid grey) from the 
envelope onto the surface of the disk (black). In the absence of 
the disk, the trajectories would extend to the midplane along the 
dotted lines. The inset in the bottom panel shows a blow-up of 
the inner 5x1 AU. 



radius. If the vertical extent of the disk is taken into account, 
the infalling material hits the disk before it can flow all the way 
to the midplane. Because part of the disk often spreads beyond 
the centrifugal radius, especially at early times (Fig. [TJ, accre- 
tion now occurs at much larger radii. This is visualised in Fig. [3] 
which shows the mass loading onto the disk at 2.0 x 10^ yr (0.23 
face) after the onset of c ollapse. T he model parameters are those 
of the default model of lDAWOo : an initial cloud core mass Mq 
of 2.5 Mq, a rotation rate Qq of 1 x 10"'^ s"' , and a sound speed 
Cs of 0.23 km s"'. The centrifugal radius at 2.0 x 10^ yr is 2.2 
AU, but the disk has already spread to 32 AU. Accretion occurs 
across the entire disk, although most mass falls in at small radii. 

If the vertical structure of the disk is ignored when calcu- 
lating the source function, as happened in the old method of 
DAW06, the infall trajectories continue along the dotted lines. 
They all intersect the midplane inside of R^. In this case, that 
means all accretion takes place inside the "annealing radius" 
(^ann, the radius corresponding to 800 K) and all dust is turned 
into crystalline form. In the new method, only 36% of the ac- 
creting material comes inside of ^ann, so the disk gains a much 
smaller amount of crystalUne dust. 

The crystalline fractions obtained with the new method are 
compared to the results from DAW06 in Fig.|4] In both cases, the 
inner part of the disk, out to a few AU, is fully crystalline. This is 
followed by a near-powerlaw decrease as crystalline material is 
mixed to larger radii. At a few tens of AU, the crystallinity levels 
off to a base value that remains roughly constant to the outer 
edge of the disk, where the curves are terminated. The increase 
in crystallinity towards large radii in the original DAW06 curves 
is an artifact from the pre-infall low-density seed disk (used for 
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reasons of numerical stability) and occurs only in those regions 
where the surface density is anyway nearly zero. This increase 
is therefore iiTelevant. Due to the use of an improved numerical 
integration scheme, it is no longer present in the new curves. 

At each of the three time steps plotted in Fig. |4] the new 
crystallinity outside of a few AU is lower than the old one. For 
example, at 10 AU and 3.1 Myr, the fraction is down from 15.1 
to 7.4%. At the outer edge of the disk, the new method produces 
crystalline fractions down to 1%. The differences between the 
old and the new method become even more pronounced when 
we take a slightly lower initial rotation rate of 3 x 10"'^ s"' (Fig. 
|5]). The crystalline fraction at 10 AU and 3.1 Myr is now down 
from 72.6 to 11.3%. 

The amount of crystallisation that takes place during the for- 
mation and evolution of the disk depends on the initial condi- 
tions. IDAW06 already discussed the effect of the rotation rate 
(Qo) of the collapsing envelope. The more rapid the rotation, the 
larger the radius at which the bulk of the accretion takes place. 
This results in less material being heated above 800 K, so the 
disk becomes less crystalline. This effect also occurs in our new 
model. Two other conditions that can easily be changed are the 
initial cloud core mass (Mq) and the effective sound speed (c,). 
In order to get a first understanding of their effects, we co mputed 
the crystalline fractions for the parameter grid from lV09l as well 
as for models with initial masses of 0.2 and 5.0 Mq. Table [T] 
lists the relevant model parameters together with the fraction of 
crystalline silicates at three different positions and three different 
times. The first of the three positions (10 AU) is representative of 
the region probed by the recent Spitzer observations, while the 
other two positions (30 and 100 AU) contain colder material that 
can be studied with the Herschel Space Observatory. Results for 
models where the disk mass exceeds ~30% of the stellar mass 
should be treated with caution, because we do not include effects 
of self-gravity when computing the viscous evolution of the disk 
(IVorobvovl2010h . 

The models with a low sound speed all have a lower crys- 
tallinity than the models with a high sound speed. A lower sound 
speed gives a lower accretion rate from the envelope to the disk, 
as well as a lower accretion rate from the disk to the star. The 
net mass gain for the disk (accretion from envelope to disk mi- 
nus accretion from disk to star) is also lower. However, this is 
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Fig. 5. As Fig.S but for Oq = 3 x lO^'^ s 



more than offset by the accretion time becoming longer. At the 
end of the accretion phase, the disk mass is therefore higher for 
models with a low sound speed than it is for models with a high 
sound speed (Table [1). The lower accretion rate in the low-Cs 
models also gives a lower stellar luminosity, so the region where 
silicates can be crystallised is smaller. These effects combine to 
give smaller fractions of crystalline material throughout the en- 
tire disk. 

Changing the initial mass of the cloud core has a more com- 
plicated effect on the crystallinity. Table [1] shows several cases 
where, for models that differ only in the initial mass, the crys- 
talline fractions increase towards larger Mq, and several cases 
where they decrease in that direction. Due to the inside-out na- 
ture of the Shu (1977) collapse, models of different mass initially 
evolve in exactly the same way. However, the accretion phase of 
the higher-mass models in our grid lasts longer (t^^c '^ Mq) than 
that of the lower-mass models. This affects the crystallinity in 
two opposing ways. First, the protostar becomes more luminous 
for the higher-mass models (D'Antona & Mazzitelli 1994), so 
the region in which crystallisation takes place is larger. Second, 
the disk grows larger, so the bulk of the accretion occurs farther 
from the star. The first effect results in higher crystalline abun- 
dances in the inner parts of the disk in the higher-mass models, 
while the second effect eventually results in lower crystalline 
abundances in the outer parts (Fig. |6]l. Depending on the exact 
initial conditions, the transition from the inner to the outer disk 
in this context may lie inwards or outwards of the 10-100 AU 
region given in Table [T] or even within that region. 

Examples of two of these three possibilities can be found in 
the series of models with Qo = 10"'"* s ' and Cj = 0.19 km 
s"' (top part of Table [TJ. Going from Mq - 0.2 to 0.5 Mq, the 
fraction of material accreting close enough to the star to be crys- 
tallised is reduced, so we find smaller crystalline fractions at 10, 
30 and 100 AU: -6% for the 0.5 Mq model versus -20% for the 
0.2 Mq model. Increasing the initial mass to 1.0 Mq, the higher 
stellar luminosity increases the crystallinity at 10 AU to 7.6%. 
However, the crystallinity at larger radii suffers from the larger 
fraction of dust that remained amorphous during the accretion 
phase because the larger disk prevented it from getting closer to 
the protostar. At 30 AU, the crystalline fraction at 1.0 Myr de- 
creases from 6.0 to 3.8% when the initial mass goes from 0.5 to 
1.0 Mq. The decrease in crystallinity at 100 AU is even larger: 
from 5.5 to 1.9%. The differences in the other three series of 
models can all be explained in similar fashion. 
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Table 1. Crystalline silicate fractions for a range of model parameters. 
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5.3. Discussion and future work 

The goal of this paper is to show how treating the disk as a mul- 
tidimensional object and correctly solvin g the pr oblem of sub- 
Keplerian accretion affect the results of iDAWOd As shown in 



Figs. El and [Sj we obtain smaller fractions of crystalline silicates 
throughout the disk. This is an improvement over the old model, 
which was noted to overpredict crystallinity compared with ob- 
servations. 

A detailed parameter study is required to judge how well our 
cuiTent model reproduces all available observations. One com- 
plicating factor in such a procedure is the observed lack of cor- 
relation between crystallinity and other systemic properties such 
as the stellar lum inosity, the accretio n rate and the masses of the 
star and the disk (I Watson et al.l2009l) . The observed absence of a 
coiTelation between two observables usually translates to a lack 
of a physical coiTe lation, bu t this i s not always the case. For ex- 
ample, |Kessler;:Silaccretal] ( |2007|) showed why the crystallinity 
is not observed to be correlated with the stellar luminosity. The 
disk around a brighter protostar is warmer throughout, so the re- 
gion from which most of the silicate emission originates lies at a 
larger distance from the star, where the crystallinity is lower At 
the same time, though, the higher temperatures mean that more 
material can be thermally annealed, so the crystallinity at all 
radii goes up. The two effects cancel each other, so the observed 
crystalline fraction is not correlated with the stellar luminosity. 
Likewise, care must be taken when interpreting other observed 
non-correlations or coiTelations. 

The best starting point for a more detailed comparison be- 
tween model and observations appears to be the observed ra- 
dial dependence of the relative abundances of specific types of 
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crystalline silicates, such as enstatite and forsterite. Our model 
can be expanded to track multiple types of silicates, each with 
their own formation temperature and mechanism. First of all, 
this ma y help in expl aining the "crystallinity paradox" identi- 
fied bv .Olofsson et alj (12009; see also Sect. 15. lb . Second, it can 
address the question whether crystalline silicates are predomi- 
nantly formed by condensation from hot gas (~1200 K), by ther- 
mal annealing at slightly lower temperatures (~800 K), or by 
shock waves outside the hot inner disk. At the moment, neither 
the observations nor the models can rule out any of these mecha- 
nisms. The crystalline fractions obtained with our model suggest 
that thermal annealing followed by radial mixing must be taking 
place and must therefore be responsible for part of the observed 
crystalline silicates. A scenario in which all crystalline mate- 
rial is formed where it is observed, according to the model of 
iBouwman et al. (2008), appears unlikely. 

In addition to tracking multiple types of silicates, it ma y be 
worthwhile to investigate different collapse scenarios. The IShul 
(Il977h collapse starts with a cloud core with an r^^ density pro- 
file, but observations of pre-stellar cores usually show an r '-^ 
density profil e instead (Alves et al. 2001; Motte & Andre 2001; 
iHarvev et alJ 120031: lAndre et alJ 120041: iKandori et al. 2005^ 
Bonnor-Ebert (BE) spheres have such a density profile (.Ebera 
119551; lBonnoij[l956l) . so they have been proposed as an alterna - 
tive starting point for collapse models ("Whit worth et al.iri996h . 
The collapse of a BE sphere results in different dens ities, ve- 
locities MidJem^CTaturesthanAoseobtain^^ 
lapse (iFoster & Cheyalied 11 9931: 'M atsumoto & Hanawal l2003t 
[Banerjee et al. 2004l; IWalch' et al. 200SJ), leading in turn to difi'er- 
ent crystalline silicate abundances. However, no analytical solu- 
tions currently exist for the collapse of a rotating BE sphere, so 
we are unable to pursue this point in any more detail. 

6. Conclusions 

This paper presents a new method of correcting for the sub- 
Keplerian velocity of envelope material accreting onto an axi- 
symmetric two-dimensional circumstellar disk. Unlike the pre- 
vious corrections of Hueso & Guillot (2005) and Visser etal. 
(l2009l) . this new method properly conserves angular momentum 
and produces infall trajectories without discontinuities. The lat- 
ter is important for tracing changes in the chemical contents and 
dust properties during the evolution of the envelope and disk. 

The disks produced with the new method are smaller than 
those produced with the old method by up to a factor of ten. 
Depending on the initial conditions, the disk masses are between 
100 and 200% of previously computed values (Sect. [3]). The new 
disks are a few degrees colder in the inner regions and a few de- 
grees warmer in the outer regions, resulting in lower abundances 
of CO ice (Sect. H)). By the time the system reaches the clas- 
sical T Tauri stage, at about 1 Myr, the global ice abundances 
still agree well with observations. Overall, there are no majo r 
changes in the gas-ice ratios compared with lVisser et al.l (l2009h . 

The disk w as tre ated as geometrica ll y fla t by 
iDullemond et alJ (l2006ah . As in iVisser et alJ (l2009h . we 
now also take into account the vertical structure when com- 
puting the infall trajectories. This results in the bulk of the 
accretion occuring at larger radii. A smaller fraction of the in- 
falling material now comes close enough to the star to be heated 
above 800 K, the temperature required for thermal annealing of 
amorphous silicates into crystalline form. Therefore, the new 
method produces crystalline abundances that are lower by a few 
per cent to more than a factor of five compared to the old model. 
We now obtain a better match with observations and we argue 



that thermal annealing followed by radial mixing is responsible 
for at least part of the crystalline silicates observed in disks. 
An expanded model, which tracks specific forms of crystalline 
silicate, is required to establish in more detail the importance of 
this and other possible crystallisation mechanisms. 
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